package stef.home.arp.labos3;

import java.math.BigDecimal;

import stef.home.arp.labos1.matrix.Matrix;
import stef.home.arp.labos1.matrix.Op;
import stef.home.arp.labos1.utils.Utils;

public class NumericIntegration {

	private Matrix A;
	private Matrix B;
	private Matrix x;
	private double T;
	private int t;
	private int printT;
	
	public NumericIntegration(Matrix A, Matrix B, Matrix x, double T, int t, int printT) {
		super();
		this.A = A;
		this.B = B;
		this.x = x;
		this.T = T;
		this.t = t;
		this.printT = printT;
	}
	
	public Matrix[] rungeKutta(){
		System.out.println("----------------------Runge-Kutta:----------------------");
		int interval = (int)(t/T);
		Matrix[] X = new Matrix[interval+1];
		Matrix x = this.x;
		X[0] = x;
		for(int i = 1; i<=interval; i++){
			Matrix m1 = Op.additionMatrix(Op.multipyMatrix(A, x), B, Op.ADD);
			Matrix m2 = Op.additionMatrix(Op.multipyMatrix(A, Op.additionMatrix(x, Op.scalarMultiplication(T/2d, m1), Op.ADD)),B,Op.ADD);
			Matrix m3 = Op.additionMatrix(Op.multipyMatrix(A, Op.additionMatrix(x, Op.scalarMultiplication(T/2d, m2), Op.ADD)),B,Op.ADD);
			Matrix m4 = Op.additionMatrix(Op.multipyMatrix(A, Op.additionMatrix(x, Op.scalarMultiplication(T, m3), Op.ADD)),B,Op.ADD);
			Matrix temp1 = Op.scalarMultiplication(2d, m2);
			Matrix temp2 = Op.scalarMultiplication(2d, m3);
			Matrix temp = Op.additionMatrix(m1, temp1, Op.ADD);
			temp = Op.additionMatrix(temp, temp2, Op.ADD);
			temp = Op.additionMatrix(temp, m4, Op.ADD);
			temp = Op.scalarMultiplication(T/6d, temp);
			x = Op.additionMatrix(x, temp,Op.ADD);
			X[i] = Utils.copyMatrix(x);
			if(i%printT == 0){
				System.out.println(i+" iteracija:");
				Utils.printMatrix(x);
			}
		}
		return X;
	}
	
	public Matrix[] trapez(){
		System.out.println("----------------------Trapez:----------------------");
		int interval = (int)(t/T);
		Matrix[] X = new Matrix[interval+1];
		Matrix x = this.x;
		X[0] = x;
		Matrix U = Utils.createIdentityMatrix(A.getNcolumns());
		Matrix temp = Op.additionMatrix(U, Op.scalarMultiplication(T/2d, A),Op.SUB);
//		Utils.printMatrix(temp);
//		Utils.printMatrix(Op.inverz(temp));
		Matrix R = Op.multipyMatrix(Op.inverz(temp), Op.additionMatrix(U, Op.scalarMultiplication(T/2d, A),Op.ADD));
		Matrix S = Op.multipyMatrix(Op.inverz(temp), Op.scalarMultiplication(T/2d, B));
		
		for(int i = 1; i<=interval; i++){
			x = Op.additionMatrix(Op.multipyMatrix(R, x),S,Op.ADD);
			X[i] = Utils.copyMatrix(x);
			if(i%printT == 0){
				System.out.println(i+" iteracija:");
				Utils.printMatrix(x);
			}
		}
		return X;
	}
	
	public static void main(String[] args) {
		Matrix A = Utils.createMatrix("A.txt");
		Matrix B = Utils.createMatrix("b.txt");
		Matrix x = Utils.createMatrix("x.txt");
		
		NumericIntegration integration = new NumericIntegration(A, B, x, 0.01d, 10, 1);
		integration.trapez();
	}
}
